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ABSTRACT 
Recent high resolution observations of the Galactic center black hole allow for direct comparison 
with accretion disk simulations. We compare two-temperature synchrotron emission models from 
three dimensional, general relativistic magnetohydrodynamic simulations to millimeter observations 
of Sgr A*. Fits to very long baseline interferometry and spectral index measurements disfavor the 
monochromatic face-on black hole shadow models from our previous work. Inclination angles < 20° 
are ruled out to 3cr. We estimate the inclination and position angles of the black hole, as well as the 
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electron temperature of the accretion flow and the accretion rate, to be i = 50°" 

Te = (5.4 ± 3.0) X 10i"K and M = bt^^ x IQ-^Mqyy''^ respectively, with 90% confidence. The black 
hole shadow is unobscured in all best fit models, and may be detected by observations on baselines 
between Chile and California, Arizona or Mexico at 1.3mm or .87mm either through direct sampling 
of the visibility amplitude or using closure phase information. Millimeter flaring behavior consistent 
with the observations is present in all viable models, and is caused by magnetic turbulence in the 
inner radii of the accretion flow. The variability at optically thin frequencies is strongly correlated 
with that in the accretion rate. The simulations provide a universal picture of the 1.3mm emission 
region as a small region near the midplane in the inner radii of the accretion flow, which is roughly 
isothermal and has v /vc '~ 1 — 20, where v^ is the critical frequency for thermal synchrotron emission. 

Subject headings: accretion, accretion disks — black hole physics — radiative transfer — relativistic 
processes — galaxy: center 



1. INTRODUCTION 

Due to its proximity, th e compact radio source at the 
Galactic center (Sgr A*; 'Balic k fc Brownl ITQTI ) is the 
most intensively studied supermassive black hole candi- 
date. Usi ng very long ba s eline interferometry (VLBI) 
at 1.3mm. iDoeleman et al.l (|2008l ) detected source struc- 
ture on event-horizon scales. Future measurements may 
detect the black hole shadow, providing the flrst direct 
evidence of an event horizon (Bardccn 1973; Falck e~e"t al.l 
[2OOOI) . 

S gr A* is moderately variable in t he radio a nd millime- 
ter dzhao eFap 120031: lEckart et al.l 120 08: Ma rrone et al.l 
l2008HLi et al.ll2009l : lYusef-Zadeh et al.M 20091 with order 
of magnitude flares in the IR and X-ray ( Bagano ff et al.l 
[200lt JGenzel et all[200l iGhez et al.l [20041 ). Its spectral 
energy distribution (SED) rises from the radio to a mil- 
limete r peak ( "submillimeter bump"). IBaganoff et al.l 
()2003t ) detected a much lower quiescent luminosity in 
the X-ray, while a quiescent infrared state has not 
been definitively detected. The submillimeter bump is 
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thought to arise from relativistic, thermal electrons in the 
innermost portion of a hot, thick, adyection-dominated 
accre tion flow (ADAF; iNaravan et al.l[l995t lYuan et aD 
|2003() . The radio spectrum can be explained either by a 
small amount of the internal energy being inject ed into 
nonth ermal electrons in the accretion flow (Yua n et al.l 
l2003f) or by a short, mild ly relativistic, optically thick jet 
(jFalcke fc: Marko9l2000[ ). The quiescent X-ray emission 
is dominated by bremsstrahlung from cooler electrons 
f urther from the bla c k hole , near the Bondi radius. 

iBroderick fc Loebl ()2006[ ) modeled the infrared flaring 
as being due to orbiting inhomogeneities ( "hot spots" ) in 
the inner radii of the accretion flow. Recently, correlated 
multiwavelength flares have in stead favored an adiabat - 
ically expanding blob model (jYusef-Zadeh et ani2009( ). 
which can explain the observed time lags between the 
infrared. X-ray and millimeter flares. The jet emission 
model can explain the time lags between flari ng events 
in the radio and millimeter (jFalcke et al.ll2009l ). 

The ADAF, jet, hotspot and expanding blob mod- 
els for Sgr A* emission provide a variety of competing 
pictures of the accretion flow. These models are either 
completely non-relativistic or modified later to be consis- 
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tent with relativity. They neglect the magnetic fields re- 
sponsible for outward angular momentum transport and 
accretion via t he m agnetorotational instabi lity (MRI, 
Velikhovl 119591 : [Chandrasckhar 1960; Balbus fc HawlevI 
1991*). None can self-consistently account for quiescent 
and flaring behavior. 

Magnetohydrodynamic (MHD) simulations can pro- 
vide a more physical description of the accretion 
flow. They have been used to model the syn- 
chrotron emission from Sgr A*, either in three spa- 
tial dimensions wi t h a pseudo-Newtoii i an potential 
JOhsuga ct al. 2005; G oldston et~all 120051: IHuang et all 
[2009: .Chan et al.. .200tf) or in two dimensions in full 
general relativity (iNoble et al.ll2007HMoscibrodzka et al.l 
120091 iHilburn et al.ll2009[ ). Non-relativistic simulations 
are especially inappropriate for modeling the millimeter 
emission, which originates in the innermost portion of 
the accretion flow where relativistic effects are strongest. 
Axisymmetric simulations cannot sustain the MRI, and 
cannot accurately model variability. 

The submillimeter bump is of particular interest. This 
is the frequency range where VLBI can both resolve 
the black hole event horizon, and where the measure- 
ments are no longer dominated by interstellar scatter- 
ing. Unlike the radio, IR, and X-ray emission the 
dominant emission mechanism and elect ron distribution 
function is known in the millimeter. iBroderick et al.l 
(j2009t) fit a non-relativistic, semi-analytic radiatively 
inefficient accretion ffo w (RIAF) m o del t o spectral 
and mm VLBI data. iDexter et all (|2009( ) fit time- 
dependent images of millimeter synchrotron emission 
from a th ree-dimensional GRMH D code to the VLBI 
data from iDoeleman et al.l ()2008[ ). In this work, we ex- 
pand that analysis in several respects. We produce im- 
ages from additional three-dimensional GRMHD simula- 
tions, including total energy co nserving simulations from 
iMcKinnev fc BlandfordI (120091) and sim ulations at low 
black hole spins from [Fragile et all ([20091). Images are fit 
to the spectral index measurements from Marronc f2006') 
as well as the VLBI data from Doclcman ct al. (2008), 
and two-temperature models of the accretion ffow are 
considered. 

2. SIMULATIONS 

Relativistic, three dimensional simulations are neces- 
sary to produce a physical, realistic description of the 
time-dependent accretion flow structure. Axisymmet- 
ric simulations cannot sustain turbulence and exaggerate 
variability. Previous studies using 2D simulations have 
time-averaged the fluid structure, but even then the sta- 
tistical steady state is not correct due to the decay of the 
MRI on the local orbital time. The ideal approach would 
be to use a single code to run a grid of 3D simulations 
over black hole spin, initial magnetic field configuration 
and initial torus location and geometry. However, given 
the formidable computational expense of such a project, 
we instead leverage a subset of existing 3D GRMHD sim- 
ulations from two different groups. The simulations and 
their parameters are listed in Table [TJ 

The two sets of simulations use different initial condi- 
tions and energy evolutio n equations. Simulations from 
[Fragile et al.l ([20071 [2009D evolve i nternal energy and use 
the in itial torus configuration from De Villiers fc Hawle^ 
([2003D . whereas those from .McKinnev fc Blandford 



([2009D conserve total energy and use a n initial torus 
place d much closer to the bla ck hole ([Gammie et aT 
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2003") . The two simulations from [McKinnev fc Blandford 



2009) differ slightly in black hole spin, as well as in 
initial magnetic field configuration. The simulation la- 
beled MBD uses a standard, small-scale dipolar loop. 
The MBQ simulation has a large-scale quadrupolar field, 
which leads to a stronger MRI, slightly higher magnetic 
pressures and a larger accretion rate. 

Total energy conserving simulations are probably 
more appropriate for modeling Sgr A*, since its in- 
ferred radiative effic ie ncy i s relatively low (e ~ .001 — 
.01). [Dexter et al.[ ([2OO90 computed a radiative ef- 
ficiency from numerical energy losses of e ^ .1 
for the 90h simulation, and a bolometric "luminos- 
ity" exceeding that of Sgr A* by a factor of a 
few. However, few energy conserving three dimen- 
sional GRMHD sin iulations have be e n run to date 
(Shafce ct al. 2008: [Noble et al.[ [20091: [Noble fc Krolik[ 
12009.: McKinnev & BlandfordI [2Q09[: [Noble et al.[[2010D . 
Only [McKinnev & Blandfordf ([2009[) simulated a thick, 
advection-dominatcd accretion flow. The simulations 
used here do not include radiative cooling, which is likely 
a good approximation for an ADAF.. Non-conservative 
simulations, where a significant amount of energy is 
lost to artificial numerical cooling, are still advection- 
dominated. 

For comparison, we have also run a set of axisym- 
metric simulations usin g the publicly a vailable HARM 
code (Gammie et al. 2()0l INoble et al] [2006 ) at black 
hole spins a/M = 0.00, 0.50, 0.75, 0.90, 0.92, and 0.94 
(G = c = 1 is used throughout this paper). They are 
not used for fitting to observations due to the problems 
with axisymmetric simulations described above. Instead, 
they give a sense of the spin dependence of total energy 
conserving simulations, which helps to break the degen- 
eracies between black hole spin and initial conditions in 
our 3D simulations. They also provide an opportunity 
for direct comparison with results from previous work 
using similar simulations. 

3. MODELING 

Relativistic synchrotron radiation is computed from 
the simulations by tracing rays backwards from an ob- 
server's camera at infinity through the accretion flow. 
Each ray constitutes a pixel of the image, which are 
spaced uniformly on a rectangular grid. The points along 
each ray are assumed to be null geodesies of the Kerr 
spacetime, and their tr ajectories are comput ed using the 
pubfic code geokerr ([Dexter fc Agoll 120091 ). The code 
computes all coordinates of the geodesies, and the ray 
tracing procedure is fully time-dependent. Tracing rays 
backwards is ideal for Sgr A*, where the images them- 
selves are observables. 

We find that an image resolution of 150x150 pixels 
(rays), with maximum spacing of .2M in impact param- 
eter at infinity is sufficient for total flux convergence to 
~ 1% relative to 300x300 or 600x600. The largest single 
pixel errors are ~ 5 — 10%, concentrated near the circu- 
lar photon orbit. We flnd ^ 400 points on each geodesic 
to be adequate for all models. These points are spaced 
evenly in 1/r, except near radial turning points, where 
geokerr uses the polar angle as the independent variable 
for better resolution. 
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TABLE 1 

TO THE Simulations 




Name 


Spin (M) 


Energy Scheme 


At 


(Used in Modeling) 




Initial Conditions 


Reference 


Oh 
50h 
90h 
MBD 


0.00 
0.50 
0.90 
0.92 
0.94 


Internal 

Internal 

Internal 

Total 

Total 




7900M {4600M) 
7900M (4600M) 
7900M (4600M) 
3500M (2000M) 
5500M (3000M) 


Gamm 


De Villiers & Hawley (2003) 

De Villiers & Hawley (2003) 

De Villiers & Hawlev (20031 

Gammie et al. (2003) 


FraEileet al. (2009) 

FraEile et al. (2009) 

Fragile et al. (2007) 

McKinney & Blandford (2009) 


MBQ 


lie et al. f2003) with quadrupolar field 


McKinney & Blandford (2009) 



The unpolarized synchrotron emission coefficient is 
computed from Leung et al (2010; submitted) by linearly 
interpolating simulation fluid variables to the points 
along each ray. For reasons discussed below, we often use 
data from the nearest time step rather than interpolating 
in the time direction. The absorption coefficient is ob- 
tained from Kirchoff's Law. The gas pressure, magnetic 
field and particle density are taken from the simulations 
and scaled to physical units with a single free parameter 
representing the total torus mass. This parameter fixes 
the mass accretion rate. Only emission from fluid within 
r = 25M is used for computing images, and only the last 
^ 60% of the simulation data is used, after transients 
from the onset of turbulence have died down (see Table 
[2]). The millimeter emission in the models considered is 
always concentrated well within both our radial bound- 
ary and the initial pressure maximum for all simulations. 

The simulations evolve the total gas pressure of elec- 
trons and ions, while only the electron temperature is 
used in calcula t ing th e emergent radiation. We f ollow 
IGoldston etHI (|2005D : IMoscibrodzka et alj ([2009) and 
generate images using two-temperature models with con- 
stant ratios of Ti/T,,, which is then another free param- 
eter. 

The models a re fit to both t h e VL BI observations 
at 1.3mm from Doele man et al.l (|2008l ) and the spec- 
tral index meas urements between .4mm and 1.3mm from 
iMarronj (|2006( ). This gives the probability of observ- 
ing the measured value s from a particular model. Fol- 
lowing |BrodericF^eF^L| (|2009[ ). this distribution is con- 
verted to the probability distribution as a function of 
our parameters given the observations using Bayes' theo- 
rem. The result is two separate probability distributions: 
Pv{Mv, it,, ^, Ti/Te, a) and Ps{Ms-, ts, Ti/T^, a), where t is 
the observer time, M is the time-averaged accretion rate, 
^ is the orientation of the black hole spin axis projected 
on the sky, Ti/T^ is the ion electron temperature ratio, 
and a is the spin value of the simulation used. The sub- 
scripts V and s refer to the VLBI and spectral fits re- 
spectively, and the simulations differ in more respects 
than just the value of the black hole spin. To produce 
parameter estimates, we average over observer time and 
marginalize over accretion rate and multiply the result- 
ing distributions. To estimate the value of a single pa- 
rameter, we then marginalize the combined probability 
distributions over the other parameters. To estimate the 
accretion rate, we average over observer time, combine 
the resulting distributions, and then marginalize over all 
other parameters. 

For each simulation, images are produced over a grid of 
the above parameters at ■4mni and 1.3m m. VLBI fits are 
done as in lDexter et al.l (j2009( ) following iBroderick et al.l 
()2009f ). Visibility amplitudes are calculated by taking 



the absolute value of the Fourier transform of the 1.3mm 
image averaged over 10 minute intervals to match the ob- 
servations. The visibility is then rotated to the desired 
sky orientation, and multiplied by an elliptic al Gaussian 
to account for interstellar sc attering as done in l Fish et al.l 
(|2009a) using the fits from lBower et al.l (|2006D . The ro- 
tated, scatter-broadened visibility amplitude is then in- 
terpolated to the baseline locations of the VLBI mea- 
surements, and fit to their measured values. 

The total fluxes at .4mm and 1.3mm are time-averaged 
over 2.5hr inte rvals to mimic the observations from 
iMarronel (|2006D . We average the observed fluxes at 
1.3mm and spectral indices between .4mm and 1.3mm 
over the four observational epochs, and fit the theoret- 
ical values to the averaged observational values. Un- 
certainties are determined by adding half the range of 
the individual measurements to the instrumental errors 
in quadrature. This procedure gives a 1.3mm flux of 
3.75±.48 Jy and a spectral index of — .18±.34, where the 
uncertainties are one sigma. Note that the 1.3mm flux 
during these observations was ^ 50% higher than durin g 
the VLBI measurements from iDoeleman et alJ ()2008fl . 
Image intensities are converted to fluxes using a black 
hole mass of 4 x IQ^Mq and a distance of 8 kpo 

Since the individual measurement errors from IMarronel 
are smaller than the flux differences between 



the measurements, the data provide information about 
source variability as well as a time-averaged value. In- 
stead of using an averaged value as described above, we 
could instead flt to each measurement separately. Our 
light curves have been analyzed both ways. Using this 
second method provides slightly tighter constraints for 
individual models, but the parameter estimates and best 
flt models discussed below are identical between the two 
methods. 

4. RESULTS 

The model of the millimeter emission from Sgr A* as 
realized in the simulations is of a compact, hot, mag- 
netized accretion flow in the vicinity of a black hole. 
The millimeter is the peak of the spectrum and the fre- 
quency range which produces the majority of the ra- 
diation. It is produced by thermal electrons near the 
midplane of the disk in the inner radii, with typical 
values of n ~ 5 X lO'^cm"^, B '^ 50G, and Tg ^ 
5 X 10^°K. These typ ical values have been found by a 
number of modelers (I Yuan etalj [ 20031 : IGoldston et alJ 
l2005 HMoscibro dzka et al."2009l. 

As found in [Dexter et al. (2009<), the GRMHD simu- 
lations provide excellent fits to the observed millimeter 
emission from Sgr A* (reduced x^ < 1.5 frequently for 
all simulations). This is still true when we incorporate 
spectral index information as described above. However, 
the inclusion of that observational constraint and of addi- 
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Fig. 1. — Normalized probability density as a function of black 
hole spin, marginalized over inclination angle and Tj /Te . Tr i angles 
denote the first three entries in Table [l] from [Fragile et al.lJI2007l . 
[2009 ). Di amonds repre s ent th e last two entries in Table [T] from 
iMcKinnev fc BlandfordI 1(20091) . 
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Fig. 4. — Normalized probability density as a function of electron 
temperature, marginalized over inclination angle and simulation 
(black hole spin). To create this plot, the grid of simulations over 
Ti/Te was binned by median electron temperature in the region of 
largest emissivity. Probability densities for models falling in each 
bin are added to make the histogram. The dotted histogram is the 
result of assigning equal probability to each model. 
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Fig. 2. — Normalized probability density as a function of ob- 
server inclination, for each simulation separately and marginalized 
over simulation (black hole spin). All curves are marginalized over 
Ti/T^. The Oh simulation isn't shown due to its negligible contri- 
bution to the overall curve, but is similar to that of 50h. 
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Fig. 3. — Normalized probability density as a function of sky 
orientation, for all simulations (solid) as well as each simulation 
individually. 
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Fig. 5. — Normalized probability density as a function of accre- 
tion rate, marginalized over inclination angle and simulation (black 
hole spin). Probability densities for models with accretion rates 
falling in each bin are added to make the histogram. The dot- 
ted histogram is the result of assigning equal probability to each 
accretion rate sampled. 

tional simulations allow us to estimate some parameters 
of the accretion flow. Single parameter probability distri- 
butions are found by marginalizing the combined spectral 
and VLBI probability distributions over the others. To 
estimate c onfidence intervals, we use the technique de- 
scribed in IBroderick et all ()2009( ). The probability dis- 
tribution is integrated from its maximum value, Pmax, 
down to a cutoff value Pmin so that the cumulative prob- 
ability enclosed between Pmin and Pmax is equal to the 
desired confidence interval. Then the estimated value, 
lower and upper bounds of the parameter x are x{pmax), 
min[x{p > Pmin)], and max[a;(p > Prmn)] respectively. 

4.1. Parameter Constraints 

We estimate the electron temperature in the millimeter 
emission region, the accretion rate onto the black hole, 
the inclination angle of the accretion flow relative to the 
observer and the projection of the black hole spin axis 
on the sky. Unfortunately, our use of existing Simula- 
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Fig. 6. — Light curves for best fit parameters (see Table [2j from each simulation at l.Smm (solid) and .4mm (dotted). The vertical 
arrow denotes the time of best fit to the VLBI measurements, while the horizontal arrows denote the averaged interval corresponding to 
the best fit to the spectral index observ ations. The a ccretion rate used for the plots is not necessarily the best fit, since that value can 
differ between the two constraints. Thc iMcKinnev ^~B landford (i200ft l simulations were run for shorter times since the torus was initially 
placed at smaller radius. 

estimate the inclination angle as i = 50°1^5o with 90% 
confidence. This interval may change with more com- 
plete sampling of the spin and initial conditions param- 
eter space; however the similarity of the curves from 
different simulations suggests that such changes will be 
minor. The fa ce-on, b l ack h ole shadow fits at i < 20° 
found in Dext er et al.l ()2009() are now ruled out to 3cr 
confidence. One reason is that the spectra l index con- 
straint, which was not used in Dex ter et al.l (|20091 . dis- 
favors low inclinations. Also, good fits to the VLBI data 
at low inclination occur infrequently in the MBD and 
MBQ simulations, which were not considered in the pre- 
vious study. The VLBI measurements also depend on 
the sky orientation, ^, the position of the black hole spin 
axis projected on the sky measured E of N. A plot of 
p(^) is shown in Fig. [31 The distribution is quite broad, 

and we estimate a value of ^ = — 23°^22° with 90% con- 
fidence. Both the estimated values and the probability 
distributions shown he re are in excellent agreement with 
iBroderick et"al] (|2009t) (c.f. their Fig. 7). 

For all simulations, the probability distribution over 
Ti/Te has a peak value between 1 and 10. The best 



tions precludes us from determining a constraint on the 
black hole spin. A plot of p{a) is shown in Fig. [TJ The 
large peak at high spin is deceptive, and is due to the dif- 
ferences between the simulations rather than black hole 
spin. While it appears we might be able to rule out low 
spin simulations, the Oh simulation is probably a special 
case. Low electron temperatures in that simulation cause 
it to have a large photosphere (r ~ lOM) even when 
Ti — Te, and therefore it fits the observations poorly. 
The axisymmetric HARM simulations, which are more 
similar to a 3D conservative simulation like MBD/MBQ, 
have larger electron temperatures at all spins and are 
optically thin for Ti = T^. So although the Oh model is 
ruled out, low spin conservative simulations should have 
larger electron temperatures and provide better fits. We 
do not consider the Oh simulation for the remainder of 
the paper due to its negligible infiuence on the overall 
probability distributions. 

All best fit models have similar probability distribu- 
tions in inclination angle. These are shown in Fig. [21 
along with the combined probability distribution after 
marginalizing over all other parameters. From this we 
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Fig. 7. — Light curve at .4mni (solid) and accretion rate at the in- 
ner boundary of the MBQ simulation (dotted) with best fit param- 
eters (see Table [211. Both quantities are scaled to their maximum 
value. 

fit Ti/Te rises witli MHD temperature, so that rouglily 
the same electron temperature fits best across all simula- 
tions. The probability distribution of emissivity averaged 
electron temperature is plotted in Fig. U) The dotted 
histogram shows the result of assigning equal probability 
density to each model. With better sampling, the dis- 
tribution would likely be much more smooth. The elec- 
tron temperature is estimated to be (5.4 ± 3.0) x 10^"K 
with 90% confidence. As discussed below, our sim- 
ulated light curves are consistent with an isothermal 
emission region. This electron temperature is consis- 
tent with electron temperatures fou nd in previous work 
jYuan ct al. 2003; Goldsto n eflll 120051 : Isharma et all 
[2007: Moscibrodzka et al. 2009( 1 

Similarly, the accretion rate can be constrained from 
the joint fits despite the different ranges found for best fit 
parameters for each simulation (see Table [2|). The proba- 
bility distribution over accretion rate is shown in Fig. [5l 
and we estimate its value to be M = 5^2^ x lO~^Af0yr~^ 
with 90% confidence. The probability distributions from 
the two observational constraints have been combined to 
produce the plot, despite the fact that the total flux is 
'^ 50% higher during the spectral index measurements. 
This procedure favors models with substantial variabil- 
ity. However, probability distributions from the two con- 
straints separately are similar to the one in Fig. [5l with 
slightly broader ranges of allowed accretion rates and es- 
timated values in agreement with that of the joint fit. 

4.2. Variability 

iDexter et al.l (|2009f) found that simulation 90h exhib- 
ited millimeter flares with 2.5hr duration, .5hr rise times 
and 40% — 50% amplitude, in agreement with observa- 
tions of Sgr A*. Light curves at 1.3mm (230 GHz) and 
.4mm (690 GHz) for best flt parameters from all viable 
simulations (excluding Oh) are shown in Fig. [5] All light 
curves show millimeter flares that are consistent with the 
observations. The spectral index rises during the flares, 
often above zero, due to stronger variability at 690 GHz 
where the flow is completely optically thin. Only 1 — 2 
fiares occur in each theoretical light curve, which is con- 
sistent with their observed frequency. 

The fiares are caused by increased magnetic field 
strength in the inner radii of the accretion flow where 



Fig. 8.— Last half of the 1.3mm light curve of MBQ with Ti/T^ = 
3, i = 10° (Fu) compared to the volume integrated synchrotron 
emissivity (ji,) and approximate volume integrated emissivities, nB 
and nB-^ where n is the electron density and B is the magnetic field 
strength. The emissivities are normalized to their maximum values, 
while the light curve is scaled so that it matches the synchrotron 
emissivity at ~ 11 hours. 
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Fig. 9. — Light curves of MBQ at i = 10 and 70 degrees for 
n = 1, 75 where n is the number of simulation time steps used to 
compute the intensity at a single observer time, n = 1 neglects the 
effects of finite light travel time through the fluid. 
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Fig. 10. — Light curves of MBQ at 10 degrees using simulation 
time steps separated by At = 2, 8, 20 and 50M, where 2M is the 
value used for other MBD and MBQ light curves, while Ai = 4M 
for the 50h and 90h simulations. 
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the synchrotron emissivity is highest. This is the case 
even in total energy conserving simulations, where the 
gas can be heated from reconnecting magnetic fields. Re- 
connection events, while present, do not seem related 
to the millimeter variability. Likewise, the variability is 
not produced by local inhomogeneities (hotspots), but 
rather by relatively global increases in the magnetic field 
strength in the inner radii. The light curves at .4mm are 
strongly correlated with variability in the accretion rate, 
as shown for simulation MBQ in Fig. [71 There are no 
robust lags found between the variability in the millime- 
ter light curves and the accretion rate, and we cannot 
definitively determine the cause and effect. 

Since only the millimeter emission is modeled here, it 
is unclear if the flares produced are consistent w ith the 
obser v ed x-ray, IR and mil li meter f lares (Ecka rt et al.l 
2008; 'Marrone et al.' [200l IDodds-Ed en et al.l ~ l200a 
Yusef-Zadeh et al. 20091). At 1.3mm, the flares are 



TABLE 2 
Best Fit Model Parameters 



weaker at higher inclination, where they are attenuated 
by higher optical depth between the observer and the 
emission region. Due to the presence of a significant pho- 
tosphere, the correlations with accretion rate variability 
are much weaker at longer wavelengths as well. 

The 1.3mm light curve shapes are closely reproduced 
by assuming the optically thin synchrotron emission 
comes from an isothermal region with a magnetic field 
satisfying v/v^ '- 1 - 20, where v^ = 6.27 x W^B{kTf 
(cgs) is the critical frequency for synchrotron emission. 
Then the emissivity is j^ ^ nB", with 1 < a < 2. The 
light curve from simulation MBD is compared to the in- 
tegrated synchrotron emissivity as well as these simpli- 
fied isothermal models in Fig. [51 A low inclination light 
curve is used to lessen the impacts of Doppler beaming 
and optical depth. The flaring behavior is captured with- 
out any temperature fluctuations, and the variations are 
caused by those in the magnetic field strength and par- 
ticle density. In models with significant photospheres, 
including most at high inclination, these simple formulas 
still work to describe the integrated synchrotron emissiv- 
ity, but tend to be more variable than the light curves 
which i nclude self-a bsorption. 

As in lDexter et al. (2009) . we compute radiative trans- 
fer in full general relativity, including accounting for time 
delays along geodesies. Including the time delays means 
that many different time steps of simulation data are 
used to compute the intensity along each geodesic. Fig. 
[SI shows a comparison of light curves with and without 
accounting for time delays. The light curves are virtually 
identical at both low and high inclinations, except for a 
systematic shift in observer time, indicating the constant 
time lag between r = 25 M, where we begin following 
geodesies, to the emission region. 

This is surp r ising especially given the results of 
INoble fc KrolikI (|2009D . They found significantly steeper 
power spectra at high frequencies when accounting for 
time delays. However, the simulation from that work 
used an artificial emissivity, lower time resolution (At = 
20il/) and only a 7r/2 wedge in (/). We've tried match- 
ing their 7r/2 wedge in </) and sampling with the same 
time resolution. In all cases, we find no systematic dif- 
ferences when accounting for light travel time through 
the accretion flow, and the errors due to neglecting time 
delays are at the 10% level (Fig. [9]), comparable to the 
overall uncertainty in the interpolation scheme or emis- 
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sivity used. The largest errors occur at high inclina- 
tion, where the light travel time to different parts of the 
emission region is largest. When using both the limited 
domain in (f) and At — 20M, the slope of the power spec- 
trum is slightly steeper when using time delays. However 
th e deviation is considera bly smaller than that reported 
bv INoble fc KrolikI (|2009f ). possibly due to the differing 
emissivities or disk scale heights. They also included 
cooling self-consistently in the simulation. 

Their use of a 7r/2 wedge in cf) is predicated on the as- 
sumption that the dominant azimutha l spatial structure 
is in modes with to ~ 4 and higher (jSchnittman et aLl 
12006( 1. We test this by expanding fluid variables in spher- 
ical harmonics for different simulations and time slices. 
In all cases, we find that the m = 1,2 power is larg er than 
at higher orders, as found previously by Henise v et al.l 
(2009) for a tilted simulation. This suggests that the use 
of restricted domains incorrectly approximates the spa- 
tial structure of a fully global simulation. In particular, 
we find that the power law index of the power spectrum 
is significantly steeper when using the full 27r, indicating 
that using a limited wedge in introduces spurious vari- 
ability on shorter timescales. Studying variability using 
simulations probably requires using the full (j) domain. 

We also caution that interpolating magnetic field vec- 
tors between time steps with insufficient time resolution 
results in systematic suppression of the total flux. This 
occurs when magnetic timescales are significantly shorter 
than the time step between simulation data dumps, caus- 
ing interpolation between uncorrelated magnetic field 
vectors. For the MBD and MBQ simulations, 2 - 8M 
time resolution is sufficient, as shown in Fig. [TOl How- 
ever, when interpolating magnetic fields between time 
steps in the 90h simulation, the fiux is systematically 
^ 10% lower with a time spacing of 4M. Our radiative 
transfer code uses time steps equal to an integer multiple 
of the simulation time step so that these errors cannot 
cause spurious features in the light curves. Regardless, 
we use data from the nearest time step rather than in- 
terpolating magnetic fields when the time steps are too 
large. 

4.3. Best Fit Models 

Table[2llists the parameters for best fit models from the 
viable simulations (excluding Oh), as determined by the 
joint spectral index and VLBI fits. Images and visibili- 
ties from these models at 1.3mm and .87mm are shown 
in Fig. [TTl at the time of best fit to the 1.3mm VLBI 
observations. All best fit images are essentially the su- 
perposition of two crescents, due to the Doppler beaming 
of gas approaching the observer and the circular photon 
orbit. The shadow of the black hole is prominent along 
the vertical axis as a local minimum in the visibility am- 
plitude. One dimensional visibility profiles, plotted along 
the line of the average location of the long baseline mea- 
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Fig. 11. — Images and visibility amplitudes at 1.3mm (first two rows) and .87mm (bottom two rows) for the best fit models. The first and 
third rows are images, while the second and fo urth are the corresponding visibility amplitudes. All are rotated to their best fit positions. 
The uv-plane locations of the baselines used in Docleinan ct al. (2008) are over-plotted on the visibilities. The panel size is I50fj,as, and 12 
GX for the visibilities. The columns are labeled by simulation, and each panel is scaled to its maximum value. At 1.3mm, the maximums 
are always ~ 2.4Jy. However, the total flux can vary between simulations at .87mm. Since the images and visibilities form a Fourier 
transform pair, in general a larger image corresponds to a smaller visibility and vice versa. 



surcments from lDoeleman et al.l()2008[ ). are shown in Fig. 
[T2l In all models, the visibility profile decreases mono- 
tonically with baseline length along this axis. Therefore, 
the black hole shadow may not be see n by future mea- 
surements using the same telescopes as iDoeleman et alJ 
(gOOS). 

The shadow is detectable from the two dimensional 
visibilities at baselines of similar length and roughly or- 
thogo nal orientation to those used by IDoeleman et al.l 
(|2008| ) in all best fit models. At .87mm, the resolution 
of the array improves, but the images are also smaller. 
Depending how optically thin the model is (how large 



and negative the spectral index is) , the image can be ex- 
tremely small as is the case f or MBQ. For a sp ectral index 
similar to that observed bv iMarroM ()2006[ ). the image 
will still be large at .87mm, and the shadow should be 
observable from the visibility. 

To demonstrate this explicitly, visibility amplitude 
as a function of baseline length is plotted in Fig. 
[T3l for the best fit models, considering baselines 
made up of combinations of telescopes in California 
(CARMA), Arizona (SMTO), Mexico (LMT) and Chile 
( APEX/ ASTE/ ALMA). The shadow is detectable for 
the best fit MBD and MBQ models on the Chile-Mexico 
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200S). In these cases, the closure phase ([Rogers et al.l 



1974f ). the sum of the phases measured along a triangle 
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Fig. 12. — One dimensional visibi lity profiles along the li ne of the 
average long baseline location frorn'Doelcman et aLIBOOSI ) . Unlike 
the best fit models in Dexter ct al. (2009), all of these profiles de- 
crease monotonically with baseline length. 
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Fig. 13. — Visibility amplitude as a function of baseline length 
for best fit models at 1.3mm. The telescopes considered are in 
Arizona (SMTO), California (CARMA), Mexico (LMT) and Chile 
(APEX/ASTE/ALMA). The black hole shadow appears as a min- 
imum in the ID profile. 



and Chile- Arizona baselines. The differences in the pre- 
dictions from the models in Fig. [T3]are almost exclusively 
caused by their different best fit sky orientations. 

Since Sgr A* varies significantly on timescales short 
compared to its mutual visibility between current and fu- 
ture VLBI telescopes, the profiles shown in Fig. [T51will 
necessarily include source variability. In addition, the 
short coherence time of the atmosphere and varying sig- 
nal/noise between telescopes in the arrays may prevent 
proper calibration of the visibilities (.Doclcman ct all 



of baselines, may be used as a non-imaging observable. 
It is independent of calibration errors, and probes the 
source structure in the uv-plane. IDoeleman et al.l ()2009f ) 
used a RIAF-t-hotspot model to compute predicted clo- 
sure phases for future VLBI experiments. We make pre- 
dictions for the best fit models discussed here, whose 
time-dependent structure comes directly out of the tur- 
bulence driven by the MRI. 

Fig. [TJ] shows the closure phase maps vs. observer 
time and source sky orientation for the current Ari- 
zona/California/Hawaii triangle and a future one con- 
sisting of Arizona/Mexico/Chile. Due to the moderate 
inclination favored by the VLBI and spectral observa- 
tions, the closure phase is relatively robust over time, 
except over a narrow range of sky orientations. At these 
positions, one of the baselines in the triangle is sampling 
the black hole shadow, which causes a sudden rotation 
of the phase when the real part of the visibility switches 
sign. The earth's rotation moves the baseline positions, 
slowly changing the sky orientations at which this effect 
is visible. For large phase rotations, the visibility ampli- 
tude is usually quite small (< .IJy) and will be difficult to 
detect. However, it could be observed as a transition in 
phase over a short timescale, where Sgr A* is detected at 
one fairly stable closure phase, goes to an undetectable 
amplitude, and reappears at a different stable closure 
phase. 

If the sky orientation of Sgr A* is ^ — 50°, this 
shadow signature would be detectable using the Ari- 
zona/Mexico/Chile triangle. The best fit images for 
MBD, MBQ and 90h in Fig. [13] are in this range. Note 
that the exact value of the phase rotation and how of- 
ten it occurs in time are mo del- dependent - in Fig. [14] 
it is usually ^ tt/2 and occurs at 1 — 3 observer times. 
Regardless of the value of the position angle, the black 
hole shadow can be detected from closure phase measure- 
ments on some triangle of baselines in this way. Closure 
phases from a larger set of baselines will constrain the 
structure of the accretion flow in a similar fashion to the 
current measurements of the visibility amplitude. For 
simplicity, the closure phases shown here are for almost 
co-linear baseline trios. Much more complex and time- 
dependent closure phase signatures are found on nearly 
triangular baseline arrangements, such as those including 
Hawaii, Chile and California/ Arizona/Mexico. 

We calculate the standard deviation of the small sam- 
ple of best fit visibilit ies, as done for a large ensemble 
of RIAF images bv Fi sh et all ([2009aD . The resuh from 
the best fit 3D GRMHD models, shown in Fig. [TS] is 
in relatively good agreement with their Figure 2. The 
best baselines for constraining the accretion flow are of 
similar to slightly shorter length as the existing measure- 
ments (^ 2000 — 4000km), but at roughly orthogonal ori- 
entation. The peak values of the variance are larger at 
.87mm. This is expected since we fit to existing mea- 
surements at 1.3mm, while at .87mm only the total flux 
is somewhat constrain ed by the spectr al index between 
.4mm and 1.3mm from iMarrond ([20061 ). 

The radio spectra from all models are shown in Fig. 
\W[ The shaded gray regions show the extent of the vari- 
ability at each frequency. Below ~ lO^^Hz, a substantial 
portion of the emission is produced outside r = 25M, 



10 



Dexter et al 





6 
5 
4 
3 
2 
1 



i 

1 


MBD 

.•X jm..; 



■ 




MBQ 



-150-100 -50 50 100 150 
f (degrees) 



-150-100 -50 50 100 150 
{ (degrees) 



-150-100 -50 50 100 150 
{ (degrees) 



" 3 
2 



•m 




6 r 

5 '■ 

"m" 4 r 

D 
O 

£• 3 : 
2 ■ 




■'^ 1 i 

5 f "^ ■: 

4^^ "-- 

2^4 ^ 
1 1 

1 r - 



-150-100-50 50 100 150 
f (degrees) 



-150-100 -50 50 100 150 
f (degrees) 



-150-100 -50 50 100 150 
f (degrees) 



-150-100 -50 50 100 150 
f (degrees) 



Fig. 14. — Closure phase as a function of observer time and sky orientation for best fit models on triangles of baselines between 
Arizona/California/Hawaii (top) and Arizona/Mexico/Chile (bottom). 



and so is not included here. In addition, we do not model 
emission from nonthermal electrons, and so do not expect 
to be able to fit the observations outside of the millimeter 
peak. The variability increases with frequency, due to the 
decreasing optical depth. The 50h simulation produces 
a large amount of variability, which may be inconsistent 
with the range observed from Sgr A*. 

In addition to VLBI observations, there are constraints 
on the time-dependent structure of the accretion flow 
from astrometry. Th e centroid wander of Sgr A* at 7mm 
was constrained by IReid et all (j2008[) to be < lOO/ias 
on timescales of a few hours. Fig. [T7] shows x and y 
centroid positions as functions of time for the best fit 
models at 345GIIz. The position wander for all mod- 
els considered here is < 30/xas at millimeter frequencies 
on all timescales, well within the observed upper limit. 
A prediction of the simulations is that the millimeter 
position wander shouldn't exceed that level, even dur- 
ing flares. This is due to the robustness of the time- 
dependent image structure at all inclinations, and a re- 
sult of the fact that the variability is produced globally 
in the inner radii, rather than by small, orbiting inhomo- 
geneities (hotspots). 

The instrument GRAVITY at the Very Large Tele- 
scope (VLT) will provide ~ 15/ias precision astrometry 
of Sgr A* in the infrared (IR). If the IR flaring is due 
to a relatively global process such as Compton scatter- 
ing from millimeter seed photons, centroid measurements 
from thermal synchrotron emission in the IR may be sim- 
ilar to those during the flares. Since the position wander 
from our models rarely exceeds the detection limit of 
GRAVITY at any frequency, it should be easy to distin- 
guish such IR flare mechanisms from nonthermal events 
involving hotspots or expandi ng blobs, which would lea d 
to larger centroid movements (jZamaninasab et al.ll201ClD . 



The millimeter emission region is in the innermost radii 
of the accretion flow, peaked in the midplane. The emis- 
sion is domina ted by sync hrotron emission from ther- 
mal electrons (|Yuan et al.l [20031 . The saturated mag- 
netic field strengths are typically ^ .01 — .1 of equiparti- 
tion, and both the particle density and magnetic pressure 
are proportional to the accretion rate. These quantities, 
then, are scaled together to produce the observed total 
flux from Sgr A*. For these reasons, the quantity that 
differs most between simulations is the electron temper- 
ature. The radial distribution of particle density and the 
electron temperature in large part determine the image 
structure and viability of the model. 

As expected for accretion disks with aligned angular 
momentum and black hole spin axes, the effective inner 
radius of the disk in all simulations moves closer to the 
black hole with increasing spinQ However, the density 
at the event horizon is usually only a factor ~ 5 lower 
than at the peak. There is no sharp inner boundary to 
the accretion flow. 

4.4.1. Optical Depth Effects 

For best fit parameter combinations, the 1.3mm photo- 
sphere is quite compact and the flow is mostly optically 
thin, especially at low inclination. In the Oh and 50h 
simulations, however, as well as for Ti/T^ = 10 in all 
models, the photosphere at 1.3mm extends well outside 
of the region of peak emissivity, at r ~ 5M. In general, 
the photosphere gets smaller at higher electron tempera- 
tures, since those models require lower accretion rates to 
match the observed flux from Sgr A*. Optical depth ef- 
fects are also more important at high inclination. This is 
because Doppler beaming increases the absorption along 
rays between the observer and the approaching gas. 

Single temperature conservative models (MBD, MBQ) 
have no noticeable photosphere. This leads to a small 



4.4. Emission region properties 



^ For the misaligned case, see IFragiH l|200g| ) . 
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image size, and low accretion rates. Since the iMarronj 
|2006) observations find that the turnover in the syn- 
chrotron spectrum occurs between 1.3mm and .4mm, 
they are inconsistent with models that are completely 
optically thin at 1.3mm. The VLBI observations also 
disfavor such models due to their small image sizes. For 
Ti/Tf. = 10, we have the opposite situation. VLBI rules 
out these models due to their large photospheres and 
hence image sizes, while they are generally also optically 
thick enough that the synchrotron turnover occurs short- 
ward of .4mm. For all simulations, the best fit Ti/T^. 
occurs when there is a photosphere inside of r ~ 5M. 

4.4.2. Comparison to RIAF Models 

There are significant differences between the simulated 
and semi-analytical RIAF accretion disks. In the RIAF 
model, the disk is assumed to have a constant scale height 
H/R — 1, whereas for all simulations considered here 
H/ R ~ .1 — .3 in the regions of peak millimeter syn- 
chrotron emissivity, as measured by fitting a Gaussian 
to the particle density at the median radius of peak syn- 
chrotron emissivity. The Gaussian fits describe the ver- 
tical particle density distribution, and their parameters 
are stable over time. The emissivity has the same vertical 
profile. 

Thinner disks have higher particle density, stronger 
magnetic fields and smaller temperatures. These effects 
are largely mitigated in the Sgr A* models, since the 
electron temperature scaling is a free parameter, and the 
magnetic pressure and density scale with accretion rate. 
The parameter /3 — Pg/Pm is independent of accretion 
rate and electron temperature. In regions of peak syn- 
chrotron emissivity, it has a median value around ~ 3— 10 
for the high spin simulations 90h, MED and MBQ. The 
lower spin, non-conservative simulations Oh and 50h have 
median /? ^ 10 — 20. In RIAF models, (3 is assumed to be 
a constant, usually around ^ 10. The trend of increas- 
ing /3 with decreasing spin is also seen in the axisymmet- 
ric HARM simulations. Although /? varies spatially in 
the simulations, the emission region is small enough that 
these variations aren't larger than an order of magnitude, 
so that constant /3 is likely a good approximation. 

We have also fit radial profiles of density, temperature 
and magnetic field strength to compare the simulations 
to the RIAF models. In general, simple power laws r~" 
describe the radial distribution of temperature and mag- 
netic field strength within the millimeter emission region. 
The power law indices for the temperature distribution 
are around a ~ 0.7, which is quite close to the RIAF 
value of 0.84. The conservative MBD and MBQ simula- 
tions have magnetic field indices of a ^ 1.1 — 1.2, while 
the Oh, 50h and 90h simulations range from a ^ 0.3 — 0.7. 
The RIAF value is 1.05. In the 50h, 90h and 2D HARM 
simulations, the magnetic field power law index increases 
with spin. 

Other than the disk thickness, the main discrepancy 
between RIAF models and GRMHD simulations is in the 
radial particle density profile. The simulations find flat 
or decreasing density profiles towards the horizon from 
a point at ^ 2rms, where rms is the marginally stable 
orbit. In RIAF models, the density increases up until the 
horizon as a fixed power law. This is especially inaccurate 
at low spin, where the Oh, 50h, and HARM simulations 
find peaks in the radial density profile at r ~ 10 — 15M, 




Fig. 15. — Standard deviations of the best fit visibility ain- 
plitudes at 1 .3mm and .87mm. The uv-plane locations of the 
IDoeleman et al., 1^2008) observations are over-plotted as trian- 
gles. Future baselines between Chile/Mexico/California/ Arizona 
are over-plotted as solid lines. 
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Fig. 16. — Millimeter spectra for the best fit models. The solid 
curves are the median values at each frequency, while the dark gray 
envelope shows the range reached during the simulations. Data 
points are from Falcko ct al. (1998) (open circles), An ct al. (20013) 
(filled diamonds) and Marronc (2006) (filled squares). The mod- 
els are fit to the Marronc (2006) data, while at lower frequencies 
the emission is dominated by emission outside of the simulation 
domain, and nontherma l emission from electrons eith er in the ac- 
cretio n flow dYuan et al.ii2003i) or in a short jet (iFalcke fc Marko3 
120001) . 



depending on the value of the spin. 

Finally, running the axisymmetric HARM simulations 
also allows a direct comparison to the sp ectral index val- 
ues found in iMoscibrodzka et al.l (|2009D for different ray 
tracing codes and interpolation schemes. Values of the 
spectral index, time-averaged as discussed in that work, 
are in good agreement with their results in most cases. 
This is e ven true when we u s e the angle- averaged emissiv- 
ity from IMahadevan et all (jl996[) . The emissivity from 
Leung et al (2010; in preparation) takes into account the 
angle 9 between the ray and the magnetic field. The 
largest differences between the emissivities are at ex- 
treme inclinations, where a factor of sin 9 in Vc causes 
the spectrum to peak at larger (smaller) frequencies at 
low (high) inclinations. Individual images and spectra 
can significantly differ between emissivities. However, 
all parameter estimates, probability distributions, light 
curves and best fit models presented here are robust to 
the choice of emissivity. 

5. DISCUSSION 
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Fig. 17. — X (solid) and Y (dotted) centroid positions vs. ob- 
server time for the best fit models a t 34 5GHz. The position wander 
is similar for all frequencies in Fig. 1161 



The submillimeter bump in Sgr A* provides a unique 
laboratory for connecting observations with theoretical 
models of black hole accretion flows. We have performed 
radiative transfer as a post-processor on the output from 
existing three dimensional GRMHD simulations and fit 
them to VLBI and spectral index observations of the 
Galactic center black hole. The simulations provide ex- 
cellent fits to the existing millimeter observations, and 
allow us to estimate the roughly constant temperature 
of the millimeter emission region and the inclination 
and sky orientation angles of the black hole. These are 
found to be T^ = (5.4 ± 3.0) x 10^" K, i = 50°15f„ and 

^ = — 23°^22° with 90% confidence. The estimates are in 
excelle nt agreement with those found by iBroderick et al.l 
(|2009| ) for RIAF models, despite significant differences 
between the dynamical models used. 

The face-on models, which p rovided excel l ent fi ts in 
the 90h simulation considered in lDexter et al.l (|2009t ). are 
ruled out to 3a due to the spectral index constraint in 
90h, and the paucity of good VLBI and spectral fits at 
low inclinations in the other simulations. For this rea- 
son, it may not be possible to detect the shadow of the 
black hole f rom visibility m easurements using the tele- 
scopes from lDoeleman et al.i (2008) . The best fit images 
from our modeling are crescents, with the shadow un- 
obscured in all viable models. The shadow may be de- 
tectable with measurements on a baseline between Chile 
and Mexico, California or Arizona either by directly sam- 
pling the visibility profile, or through observations of the 
closure phase. These same baselines are also the most 
promising for further constraining black hole parameters, 
and testing accretion flow models. Finally, our accretion 
rate constraint of M = 5^2^ x lO~^M0yr~^ is consistent 
with but independent of estimates from the observe d lin- 
ear p olarization and Farad ay ro tation measure s (lAgoll 
[2OO0I: [Quataert fc Gruzi"nwl[2000l: [Marrone et al.]|2007[ ). 

Millimeter light curves from all best fit models show 
flaring events with ~ 50% amplitudes and few hour du- 
rations at a frequency of 1 — 2 per day, consistent with 



the millimeter flaring behavior of Sgr A*. These flares 
are caused by rises in magnetic field strength near the 
midplane of the inner radii of the accretion flow, due to 
magnetic turbulence driven by the MRI. They do not ap- 
pear to be linked to heating from magnetic reconnection, 
and are accurately reproduced by a simple model assum- 
ing an isothermal emission region with v/vc ~ 20. The 
variability in the light curves for all simulations at .4mm 
is strongly correlated with that in the accretion rate onto 
the black hole. However, there is no clear evidence for a 
characteristic lag between the two. 

Due to the uncertainties in their emission mechanisms 
and electron distribution functions, we do not model the 
radio, IR or X-ray portions of the spectrum, and so are 
unable to test for the presence of correlated multiwave- 
length flares. iMaitra et al.l ()2009[ ) showed that the ob- 
served time lags between the radio and millimeter flares 
may be reproduced by a jet model, where perturbations 
in the accretion rate near the black hole expand as they 
flow outward. Since the millimeter variability is strongly 
correlated with the accretion rate, that mechanism could 
explain the correlations between radio and millimeter 
flares even if the millimeter emission is produced in an 
accretion disk rather than an outflow. 

The sample of simulations used does not span the pa- 
rameter space of black hole spin, initial torus geome- 
try and initial magnetic field configuration. This in- 
completeness may affect our parameter estimates. How- 
ever, the universality of the images, light curves, prob- 
ability distributions and emission region characteristics 
from all viable simulations considered here suggests that 
the changes will probably be minor. The disk thickness 
is especially important, since the emission region scale 
height mimics that of the disk. A thicker disk will also 
have larger MHD temperatures, and thus larger values 
of Ti/Te will be necessary to fit the observations. If the 
disk thickness can be constrained by VLBI observations, 
a constraint on Ti/T,, follows. This parameter would 
constrain the strength of ion-electron coupling in colli- 
sionless plasmas. Simulations with misaligned an gular 
momentum and black hole spin axes (iFragile et al.ir2007l : 
iFragile fc Blaes 20081 IFragile et al.ll20090 are also of in- 
terest, since they produce standing shocks which lead to 
considerably different image morphologies. Parameters 
such as the inclination angle or spin may be degener- 
ate with the tilt angle. Fits from these models will be 
c onsidered in a subseque nt publication. 

IBroderick et al.l (|2009l ) also estimated the parameters 
of Sgr A* from VLBI and spectral constraints by fitting 
semi-analytic RIAF models. The structure of p{i) and 
p(^) in their Figure 7 (middle and right panels) are in 
excellent agreement with ours in Figs. [2] & [3l Further, 
we find that future observations on baselines of similar 
length but orthogonal orientation to the existing Hawaii- 
Arizona baseline have the best chance of discriminating 
between o ur best fit models. This conclusion is in agree- 
ment with iFish et al.l ()2009aD . who point out that ideal 
baselines for such measurements would be between Chile 
and one or more of California, Arizona or Mexico. A 
simultaneous measurement of the 345 or 690 GHz flux 
during future 230 GHz VLBI observations would provide 
a signiflcant additional constraint. 

A grid of time-averaged images and spectra over black 
hole spin, inclination and Ti/T^ from axisymmetric, to- 
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tal energy conserving GRMHD si mulations were com- 
pared to millimeter observations by iMoscibrodzka et alj 
(|2009D . Their finding that Ti/Te = 1 models are incon- 
sistent with observations agrees with our results here for 
MBD and MBQ. As we have shown, for simulations with 
lower MHD temperatures, smaller values of Ti/Tf. will 
be necessary to reach the required electron temperature 
in the millimeter emission region. This is why the non- 
conservative simulations considered here (Oh, 50, 90h) 
are all best fit by small values of Ti/T^- It also suggests 
that such models are only appropriate for modeling Sgr 
A* when they produce large enough MHD temperatures 
that the electron temperature lies within the range es- 
timated here. At low spin, this seems to require either 
a hotter initial condition or a total energy conserving 
simulation. 

In all but the "best-bet" model from 
IMoscibrodzka et all (|2009( ). T,/Te = 10 is required 
for consistency with multiwavelength spectral obser- 
vations. Millimeter spectral index results from our 
own HARM simulations modeled in the same fashion 
are for the most part in good agreement with theirs, 
indicating that the results of these studies are probably 
robust to the ray tracing code and interpolation scheme 
employed. However, Ti/T^. = 10 models for the various 
HARM models, MBD and MBQ are all too large to fit 
the VLBI data due to their extended photospheres. In 
addition, we find good fits at Ti/T^ = 3 for MBD and 
MBQ at moderate to edge-on inclinations. We have not 
computed IR or X-ray emission, and therefore do not fit 



to the quiescent X-ray luminosity from iBaganofF et al.l 
(|2003D . or to the upper limits to the IR emis s ion. T hese 
limits are violated in IMoscibrodzka et al.l (|2009D for 
nearly edge-on inclinations and at high spin. Including 
them could disfavor high inclinations, but is unlikely to 
change our best fit models. 

Additional constraints may come from the degree 
and orientat i on of lin ear polarizat i on firs t detected by 
lAitken et afl ()2000f ). iHuang et all ()2009[) modeled the 
polarization from a pseudo-Newtonian model. It is im- 
portant to use relativistic simulations for such studies, 
since the degree of polarization is strongest close to 
the black hole. Future VLBI observations may make 
use of the closure phase between triangles of base- 
lines to probe changing structures in the accretion fiow 
(Doel eman et al.l 120091 ). iFish et al.l ([2009b) extended 
this work to include polarization. Both of these papers 
used a RIAF-f-hotspot model. We have computed the 
first closure phase signatures from more physically real- 
istic GRMHD simulations, and find that the signature of 
the black hole shadow is present in the predicted closure 
phases. Detailed comparison of observed and predicted 
closure phases will be able to constrain the accretion flow 
in the same fashion as direct sampling of the visibility 
amplitude. 
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